Chapter 2 Introduction
In the workshop, we’ll review:
- What is Spatial Data?
- What is the
sfframework for R?
## Warning: package 'sf' was built under R version 3.6.2
## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
To delve in further, let’s see some spatial data in action.
2.1 Load Spatial Data
## Reading layer `geo_export_aae47441-adab-4aca-8cb0-2e0c0114096e' from data source
## `/Users/maryniakolak/code/Intro2RSpatialMed/data/geo_export_aae47441-adab-4aca-8cb0-2e0c0114096e.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 801 features and 9 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -87.94025 ymin: 41.64429 xmax: -87.52366 ymax: 42.02392
## CRS: 4326
2.2 Non-Spatial & Spatial Views
## Simple feature collection with 6 features and 9 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -87.68822 ymin: 41.72902 xmax: -87.62394 ymax: 41.87455
## CRS: 4326
## commarea commarea_n countyfp10 geoid10 name10 namelsad10 notes
## 1 44 44 031 17031842400 8424 Census Tract 8424 <NA>
## 2 59 59 031 17031840300 8403 Census Tract 8403 <NA>
## 3 34 34 031 17031841100 8411 Census Tract 8411 <NA>
## 4 31 31 031 17031841200 8412 Census Tract 8412 <NA>
## 5 32 32 031 17031839000 8390 Census Tract 8390 <NA>
## 6 28 28 031 17031838200 8382 Census Tract 8382 <NA>
## statefp10 tractce10 geometry
## 1 17 842400 POLYGON ((-87.62405 41.7302...
## 2 17 840300 POLYGON ((-87.68608 41.8229...
## 3 17 841100 POLYGON ((-87.62935 41.8528...
## 4 17 841200 POLYGON ((-87.68813 41.8556...
## 5 17 839000 POLYGON ((-87.63312 41.8744...
## 6 17 838200 POLYGON ((-87.66782 41.8741...

2.3 Spatial Data Structure
## Classes 'sf' and 'data.frame': 801 obs. of 10 variables:
## $ commarea : Factor w/ 77 levels "1","10","11",..: 39 55 28 25 26 21 62 49 74 75 ...
## $ commarea_n: num 44 59 34 31 32 28 65 53 76 77 ...
## $ countyfp10: Factor w/ 1 level "031": 1 1 1 1 1 1 1 1 1 1 ...
## $ geoid10 : Factor w/ 801 levels "17031010100",..: 785 767 772 773 756 751 584 513 684 34 ...
## $ name10 : Factor w/ 801 levels "1001","1002",..: 782 764 769 770 753 748 545 443 663 266 ...
## $ namelsad10: Factor w/ 801 levels "Census Tract 1001",..: 782 764 769 770 753 748 545 443 663 266 ...
## $ notes : Factor w/ 7 levels "Half in CA 64 (Midway Airport)",..: NA NA NA NA NA NA NA NA NA NA ...
## $ statefp10 : Factor w/ 1 level "17": 1 1 1 1 1 1 1 1 1 1 ...
## $ tractce10 : Factor w/ 801 levels "010100","010201",..: 785 767 772 773 756 751 584 513 684 34 ...
## $ geometry :sfc_POLYGON of length 801; first list element: List of 1
## ..$ : num [1:243, 1:2] -87.6 -87.6 -87.6 -87.6 -87.6 ...
## ..- attr(*, "class")= chr "XY" "POLYGON" "sfg"
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA
## ..- attr(*, "names")= chr "commarea" "commarea_n" "countyfp10" "geoid10" ...
## Coordinate Reference System:
## User input: 4326
## wkt:
## GEOGCS["WGS84(DD)",
## DATUM["WGS84",
## SPHEROID["WGS84",6378137.0,298.257223563]],
## PRIMEM["Greenwich",0.0],
## UNIT["degree",0.017453292519943295],
## AXIS["Geodetic longitude",EAST],
## AXIS["Geodetic latitude",NORTH]]
2.4 Exploring Coordinate Reference Systems
Chi_tracts.moll <- st_transform(Chi_tracts, crs = "+proj=moll")
plot(st_geometry(Chi_tracts.moll), border = "gray", lwd = 2, main = "Mollweide", sub="preserves areas")
## Coordinate Reference System:
## User input: +proj=moll
## wkt:
## PROJCS["unnamed",
## GEOGCS["WGS 84",
## DATUM["unknown",
## SPHEROID["WGS84",6378137,298.257223563]],
## PRIMEM["Greenwich",0],
## UNIT["degree",0.0174532925199433]],
## PROJECTION["Mollweide"],
## PARAMETER["central_meridian",0],
## PARAMETER["false_easting",0],
## PARAMETER["false_northing",0]]
Chi_tracts.54019 = st_transform(Chi_tracts, 54019)
plot(st_geometry(Chi_tracts.54019), border = "gray", lwd = 2, main = "Winkel", sub="minimal distortion")
## Coordinate Reference System:
## User input: EPSG:54019
## wkt:
## PROJCS["World_Winkel_II",
## GEOGCS["GCS_WGS_1984",
## DATUM["WGS_1984",
## SPHEROID["WGS_84",6378137.0,298.257223563]],
## PRIMEM["Greenwich",0.0],
## UNIT["Degree",0.0174532925199433]],
## PROJECTION["Winkel_II"],
## PARAMETER["False_Easting",0.0],
## PARAMETER["False_Northing",0.0],
## PARAMETER["Central_Meridian",0.0],
## PARAMETER["Standard_Parallel_1",50.45977625218981],
## UNIT["Meter",1.0],
## AUTHORITY["Esri","54019"]]
## Coordinate Reference System:
## User input: EPSG:3435
## wkt:
## PROJCS["NAD83 / Illinois East (ftUS)",
## GEOGCS["NAD83",
## DATUM["North_American_Datum_1983",
## SPHEROID["GRS 1980",6378137,298.257222101,
## AUTHORITY["EPSG","7019"]],
## TOWGS84[0,0,0,0,0,0,0],
## AUTHORITY["EPSG","6269"]],
## PRIMEM["Greenwich",0,
## AUTHORITY["EPSG","8901"]],
## UNIT["degree",0.0174532925199433,
## AUTHORITY["EPSG","9122"]],
## AUTHORITY["EPSG","4269"]],
## PROJECTION["Transverse_Mercator"],
## PARAMETER["latitude_of_origin",36.66666666666666],
## PARAMETER["central_meridian",-88.33333333333333],
## PARAMETER["scale_factor",0.999975],
## PARAMETER["false_easting",984250.0000000002],
## PARAMETER["false_northing",0],
## UNIT["US survey foot",0.3048006096012192,
## AUTHORITY["EPSG","9003"]],
## AXIS["X",EAST],
## AXIS["Y",NORTH],
## AUTHORITY["EPSG","3435"]]

Chi_tracts.Hawaii = st_transform(Chi_tracts, 102114)
plot(st_geometry(Chi_tracts.Hawaii), border = "gray", lwd = 2, main = "Old Hawaiian UTM Zone 4N", sub="wrong projection!")
2.5 Refine Basic Map
## Warning: replacing previous import 'sf::st_make_valid' by
## 'lwgeom::st_make_valid' when loading 'tmap'

tm_shape(Chi_tracts) + tm_fill(col = "gray90") + tm_borders(alpha=0.2, col = "gray10") +
tm_scale_bar(position = ("left"), lwd = 0.8) +
tm_layout(frame = F)
Check out https://rdrr.io/cran/tmap/man/tm_polygons.html for more ideas.
2.6 Arrange multiple maps
tracts.4326 <- tm_shape(Chi_tracts) + tm_fill(col = "gray90") +
tm_layout(frame = F, title = "EPSG 4326")
tracts.54019 <- tm_shape(Chi_tracts.54019) + tm_fill(col = "gray90") + tm_layout(frame = F, title = "EPSG 54019")
tmap_arrange(tracts.4326, tracts.54019)
2.7 Interactive Mode
## tmap mode set to interactive viewing
tm_shape(Chi_tracts) + tm_fill(col = "gray90") + tm_borders(alpha=0.2, col = "gray10") +
tm_scale_bar(position = ("left"), lwd = 0.8) +
tm_layout(frame = F)## Warning in `$.crs`(gm$shape.master_crs, "proj4string"): CRS uses proj4string,
## which is deprecated.
## Warning in `$.crs`(crs, "proj4string"): CRS uses proj4string, which is
## deprecated.
## Warning in `$.crs`(gm$shape.master_crs, "proj4string"): CRS uses proj4string,
## which is deprecated.
## Warning in `$.crs`(crs, "proj4string"): CRS uses proj4string, which is
## deprecated.
## tmap mode set to plotting